home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / log.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  7.0 KB  |  269 lines
  1. /* specfunc/log.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_log.h>
  26.  
  27. #include "error.h"
  28.  
  29. #include "chebyshev.h"
  30. #include "cheb_eval.c"
  31.  
  32. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33.  
  34. /* Chebyshev expansion for log(1 + x(t))/x(t)
  35.  *
  36.  * x(t) = (4t-1)/(2(4-t))
  37.  * t(x) = (8x+1)/(2(x+2))
  38.  * -1/2 < x < 1/2
  39.  * -1 < t < 1
  40.  */
  41. static double lopx_data[21] = {
  42.   2.16647910664395270521272590407,
  43.  -0.28565398551049742084877469679,
  44.   0.01517767255690553732382488171,
  45.  -0.00200215904941415466274422081,
  46.   0.00019211375164056698287947962,
  47.  -0.00002553258886105542567601400,
  48.   2.9004512660400621301999384544e-06,
  49.  -3.8873813517057343800270917900e-07,
  50.   4.7743678729400456026672697926e-08,
  51.  -6.4501969776090319441714445454e-09,
  52.   8.2751976628812389601561347296e-10,
  53.  -1.1260499376492049411710290413e-10,
  54.   1.4844576692270934446023686322e-11,
  55.  -2.0328515972462118942821556033e-12,
  56.   2.7291231220549214896095654769e-13,
  57.  -3.7581977830387938294437434651e-14,
  58.   5.1107345870861673561462339876e-15,
  59.  -7.0722150011433276578323272272e-16,
  60.   9.7089758328248469219003866867e-17,
  61.  -1.3492637457521938883731579510e-17,
  62.   1.8657327910677296608121390705e-18
  63. };
  64. static cheb_series lopx_cs = {
  65.   lopx_data,
  66.   20,
  67.   -1, 1,
  68.   10
  69. };
  70.  
  71. /* Chebyshev expansion for (log(1 + x(t)) - x(t))/x(t)^2
  72.  *
  73.  * x(t) = (4t-1)/(2(4-t))
  74.  * t(x) = (8x+1)/(2(x+2))
  75.  * -1/2 < x < 1/2
  76.  * -1 < t < 1
  77.  */
  78. static double lopxmx_data[20] = {
  79.  -1.12100231323744103373737274541,
  80.   0.19553462773379386241549597019,
  81.  -0.01467470453808083971825344956,
  82.   0.00166678250474365477643629067,
  83.  -0.00018543356147700369785746902,
  84.   0.00002280154021771635036301071,
  85.  -2.8031253116633521699214134172e-06,
  86.   3.5936568872522162983669541401e-07,
  87.  -4.6241857041062060284381167925e-08,
  88.   6.0822637459403991012451054971e-09,
  89.  -8.0339824424815790302621320732e-10,
  90.   1.0751718277499375044851551587e-10,
  91.  -1.4445310914224613448759230882e-11,
  92.   1.9573912180610336168921438426e-12,
  93.  -2.6614436796793061741564104510e-13,
  94.   3.6402634315269586532158344584e-14,
  95.  -4.9937495922755006545809120531e-15,
  96.   6.8802890218846809524646902703e-16,
  97.  -9.5034129794804273611403251480e-17,
  98.   1.3170135013050997157326965813e-17
  99. };
  100. static cheb_series lopxmx_cs = {
  101.   lopxmx_data,
  102.   19,
  103.   -1, 1,
  104.   9
  105. };
  106.  
  107.  
  108. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  109.  
  110. #ifndef HIDE_INLINE_STATIC
  111. int
  112. gsl_sf_log_e(const double x, gsl_sf_result * result)
  113. {
  114.   /* CHECK_POINTER(result) */
  115.  
  116.   if(x <= 0.0) {
  117.     DOMAIN_ERROR(result);
  118.   }
  119.   else {
  120.     result->val = log(x);
  121.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  122.     return GSL_SUCCESS;
  123.   }
  124. }
  125.  
  126.  
  127. int
  128. gsl_sf_log_abs_e(const double x, gsl_sf_result * result)
  129. {
  130.   /* CHECK_POINTER(result) */
  131.  
  132.   if(x == 0.0) {
  133.     DOMAIN_ERROR(result);
  134.   }
  135.   else {
  136.     result->val = log(fabs(x));
  137.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  138.     return GSL_SUCCESS;
  139.   }
  140. }
  141. #endif
  142.  
  143. int
  144. gsl_sf_complex_log_e(const double zr, const double zi, gsl_sf_result * lnr, gsl_sf_result * theta)
  145. {
  146.   /* CHECK_POINTER(lnr) */
  147.   /* CHECK_POINTER(theta) */
  148.  
  149.   if(zr != 0.0 || zi != 0.0) {
  150.     const double ax = fabs(zr);
  151.     const double ay = fabs(zi);
  152.     const double min = GSL_MIN(ax, ay);
  153.     const double max = GSL_MAX(ax, ay);
  154.     lnr->val = log(max) + 0.5 * log(1.0 + (min/max)*(min/max));
  155.     lnr->err = 2.0 * GSL_DBL_EPSILON * fabs(lnr->val);
  156.     theta->val = atan2(zi, zr);
  157.     theta->err = GSL_DBL_EPSILON * fabs(lnr->val);
  158.     return GSL_SUCCESS;
  159.   }
  160.   else {
  161.     DOMAIN_ERROR_2(lnr, theta);
  162.   }
  163. }
  164.  
  165.  
  166. int
  167. gsl_sf_log_1plusx_e(const double x, gsl_sf_result * result)
  168. {
  169.   /* CHECK_POINTER(result) */
  170.  
  171.   if(x <= -1.0) {
  172.     DOMAIN_ERROR(result);
  173.   }
  174.   else if(fabs(x) < GSL_ROOT6_DBL_EPSILON) {
  175.     const double c1 = -0.5;
  176.     const double c2 =  1.0/3.0;
  177.     const double c3 = -1.0/4.0;
  178.     const double c4 =  1.0/5.0;
  179.     const double c5 = -1.0/6.0;
  180.     const double c6 =  1.0/7.0;
  181.     const double c7 = -1.0/8.0;
  182.     const double c8 =  1.0/9.0;
  183.     const double c9 = -1.0/10.0;
  184.     const double t  =  c5 + x*(c6 + x*(c7 + x*(c8 + x*c9)));
  185.     result->val = x * (1.0 + x*(c1 + x*(c2 + x*(c3 + x*(c4 + x*t)))));
  186.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  187.     return GSL_SUCCESS;
  188.   }
  189.   else if(fabs(x) < 0.5) {
  190.     double t = 0.5*(8.0*x + 1.0)/(x+2.0);
  191.     gsl_sf_result c;
  192.     cheb_eval_e(&lopx_cs, t, &c);
  193.     result->val = x * c.val;
  194.     result->err = fabs(x * c.err);
  195.     return GSL_SUCCESS;
  196.   }
  197.   else {
  198.     result->val = log(1.0 + x);
  199.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  200.     return GSL_SUCCESS;
  201.   }
  202. }
  203.  
  204.  
  205. int
  206. gsl_sf_log_1plusx_mx_e(const double x, gsl_sf_result * result)
  207. {
  208.   /* CHECK_POINTER(result) */
  209.  
  210.   if(x <= -1.0) {
  211.     DOMAIN_ERROR(result);
  212.   }
  213.   else if(fabs(x) < GSL_ROOT5_DBL_EPSILON) {
  214.     const double c1 = -0.5;
  215.     const double c2 =  1.0/3.0;
  216.     const double c3 = -1.0/4.0;
  217.     const double c4 =  1.0/5.0;
  218.     const double c5 = -1.0/6.0;
  219.     const double c6 =  1.0/7.0;
  220.     const double c7 = -1.0/8.0;
  221.     const double c8 =  1.0/9.0;
  222.     const double c9 = -1.0/10.0;
  223.     const double t  =  c5 + x*(c6 + x*(c7 + x*(c8 + x*c9)));
  224.     result->val = x*x * (c1 + x*(c2 + x*(c3 + x*(c4 + x*t))));
  225.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  226.     return GSL_SUCCESS;
  227.   }
  228.   else if(fabs(x) < 0.5) {
  229.     double t = 0.5*(8.0*x + 1.0)/(x+2.0);
  230.     gsl_sf_result c;
  231.     cheb_eval_e(&lopxmx_cs, t, &c);
  232.     result->val = x*x * c.val;
  233.     result->err = x*x * c.err;
  234.     return GSL_SUCCESS;
  235.   }
  236.   else {
  237.     const double lterm = log(1.0 + x);
  238.     result->val = lterm - x;
  239.     result->err = GSL_DBL_EPSILON * (fabs(lterm) + fabs(x));
  240.     return GSL_SUCCESS;
  241.   }
  242. }
  243.  
  244.  
  245.  
  246. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  247.  
  248. #include "eval.h"
  249.  
  250. double gsl_sf_log(const double x)
  251. {
  252.   EVAL_RESULT(gsl_sf_log_e(x, &result));
  253. }
  254.  
  255. double gsl_sf_log_abs(const double x)
  256. {
  257.   EVAL_RESULT(gsl_sf_log_abs_e(x, &result));
  258. }
  259.  
  260. double gsl_sf_log_1plusx(const double x)
  261. {
  262.   EVAL_RESULT(gsl_sf_log_1plusx_e(x, &result));
  263. }
  264.  
  265. double gsl_sf_log_1plusx_mx(const double x)
  266. {
  267.   EVAL_RESULT(gsl_sf_log_1plusx_mx_e(x, &result));
  268. }
  269.